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Abstract. Numerical simulations of the multi-phase interstellar medium have been 
carried out, using a 3D, nonlinear, magnetohydrodynamic, shearing-box model, with 
random motions driven by supernova explosions. These calculations incorporate the 
effects of magnetic fields and rotation in 3D; these play important dynamical roles 
in the galaxy, but are neglected in many other simulations. The supernovae driving 
the motions are not arbitrarily imposed, but occur where gas accumulates into cold, 
dense clouds; their implementation uses a physically motivated model for the evo- 
lution of such clouds. The process is self-regulating, and produces mean supernova 
rates as part of the solution. Simulations with differing mean density show a power 
law relation between the supernova rate and density, with exponent 1.7; this value is 
within the range suggested from observations (taking star formation rate as a proxy 
for supernova rate). The global structure of the supernova driven medium is strongly 
affected by the presence of magnetic fields; e.g. for one solution the filling factor of 
hot gas is found to vary from 0.19 (with no field) to 0.12 (with initial mid-plane 
field B =6nG). 
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1. Introduction 

The dynamical state of the interstellar medium (ISM) remains rela- 
tively poorly understood, especially in the presence of magnetic fields. 
Whilst the complexity of the topic still precludes a single all-encompass- 
ing model, there are many aspects of the problem that can now be 
addressed by well-designed numerical simulations. Recent 3D models 
of relevance include those of Korpi et al. (1999a, b), and of de Avillez 
(2000), de Avillez & Mac Low (2002). These models focus on small 
Cartesian boxes, with typical side-lengths of order 1 kpc, rather that 
addressing the full galactic disc. For similar reasons of computational 
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resolution, none attempt to model all of the phases of gas present, from 
the hottest ionized plasma (with density p ~ 10~ 27 g/ cm 3 , temperature 
T ~ 10 6 K) to the cold molecular clouds (p > 1CT 21 g/cm 3 , T < 10 K). 
And each of the models neglects or simplifies some aspects of the physics 
of the problem; at this early stage of study, that is to be expected. 



2. Method 

The current model — a development of the earlier work of Korpi et 
al. ( 1999a, b), already partially reported in Shukurov et al. (2003) - 
focuses in particular on the dynamics of the hot ISM gas, under the 
interaction of supernova (SN) driving with classical magnetohydrody- 
namic (MHD) forces. We model a relatively small box (of horizontal 
cross-section (0.25 kpc) 2 ), symmetric about the galactic mid-plane (z = 
0) ; this is nevertheless large enough to follow the evolution of turbulent 
ISM structures on observable scales. We neglect the very coldest gas, 
focusing on the warm, hot and ionized phases (T > 500 K). This choice 
retains enough cold gas to fit a major fraction of the mass into a small 
fraction of the volume — allowing interesting 3D structures to form 
- but allows us to run with a relatively coarse resolution of 4pc. 
Consistent with this restriction, we do not incorporate self-gravity into 
our solution, instead assuming the fixed vertical gravity profile from 
Ferriere (1998). Beyond the fine details, the numerical model is essen- 
tially as described in Brandenburg et al. (1995). The treatment of the 
thermodynamics, adapted for the ISM problem, is largely as in Korpi 
et al. (1999a,b); this approach is indebted to earlier work in 2D (e.g. 
Passot, Vasquez-Semadeni k, Pouquet, 1995; Gazol-Patiho & Passot, 
1999). 

The radiative cooling function, A, is that given by Rosen, Bregman 
& Norman (1993) (after that of Chiang & Bregman, 1998), truncated 
at 500 K to avoid the gas evolving towards the colder phases we wish to 
neglect. Heating by UV absorption, Tuv> is imposed with a net heating 
rate of 0.05 erg g -1 s^ 1 for temperatures T < 10, 000 K (e.g. Wolfire et 
al., 1995). This is not a major energy input, but again helps to prevent 
gas from cooling and contracting into phases we do not aim to resolve. 
We assume a fixed fraction of Helium of 8% (by number), and a fixed 
ionization fraction of 0.7, in implementing these two parameterisations. 

The principal driving is the supernova heating Tsn; which is im- 
plemented by the occasional, instantaneous, input of .Esn = 10 51 erg. 
This energy is locally distributed with a normalised profile of form 
exp(— (r/<i) 6 ), at distance r about the SN location; this profile, with the 
width taken as d = 16 pc, permits adequate numerical resolution. The 
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internal energy then evolves naturally; some remains in thermal form or 
is lost to radiative cooling, some is transformed to kinetic and magnetic 
energy by normal gas dynamic and MHD processes. We find this scheme 
works well, with one refinement. In gas above a certain density, injecting 
Esn in this way distributes the energy over a significant total mass; 
the central gas is then insufficiently heated, compared with the point- 
source energy injection we are trying to model. The resultant state 
undergoes rapid radiative cooling, losing too much energy this way, 
and converting an unrealistically small amount to kinetic energy by 
pressure-driven expansion. To avoid this problem, we model the initial 
phase of SN expansion by creating a cavity of width d of low density 
gas, and moving the original mass to a symmetric shell beyond this. 
Tests on the subsequent evolution show that Sedov-type expansion is 
well replicated, and both energy and mass are well conserved. With the 
SN heating implemented in this way, about 15% of the injected energy 
is typically transformed to kinetic energy. 

We determine SN frequency and location using the following scheme 
from Gudiksen (1999). (Although again, the method has antecedents: 
e.g. Bania & Lyon, 1980; Passot et al., 1995.) We identify regions with 
p > 1CT 24 g/ cm 3 , T < 4000 K as clouds, and only permit SN II within 
such clouds. The mass of gas in such regions, M c , is monitored, and 
the instantaneous SN II frequency, ^snii, is taken as 

M c X*Xst<s 

^SNII = T7 • 

Here X* = 0.02 is the gas mass fraction converted into stars (cf. Wilson 
& Matthews, 1995); X^tss = 0.1 is the fraction of the stellar mass in SN 
progenitors, obtained from a suitable initial mass function (e.g. Padoan 
et al, 1997); M S n = 10 M is the SN progenitor mass; and r c = 20Myr 
is the appropriate gas recycling time (or 'cloud lifetime'). We emphasise 
that this formula, derived on physical grounds, results in a SN driving 
that is self-regulating, with feedback from the evolving 3D structure; 
the simulations evolve towards some equilibrium SN II rate, which be- 
comes a testable output rather than a fixed input parameter. (Although 
some calibration was performed to obtain an appropriate value for r c ; 
see below.) The location of the SNe is chosen from those positions 
satisfying the criteria above, with the probability at any position being 
proportional to the local density; SN II thus end up concentrated in the 
densest regions, as is realistic. Our model also retains an appropriate 
(small) number of Type I SNe, implemented with fixed temporal and 
spatial probabilities, as in Korpi et al. (1999a); these do not play a 
large role in the solution, however. 



sarsonetal.tex; 2/02/2008; 5:15; p. 3 



4 



G. R. Sarson et al. 



Table I. Three models with varying initial mid-plane density, po. The filling 
factors are for gas with T > 10 5 K, within \z\ < 0.2 kpc. The initial magnetic 
field at the mid-plane is Bo = 6 pG in all cases. 
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Figure 1. Time-evolution of the SN II rate (left) for three simulations of differing 
density: the Arm (solid) , Average (dashed) and Interarm (dotted) models described 
in Table I. The power law fit of the mean SN II rates with density (right). 



3. Results and Discussion 



We have run simulations with three different mean densities, chosen to 
model the ISM in regions within and between spiral arms, as well as 
for an average value; these are identified as the Arm (high p), Inter- 
arm (low p) and Average cases, and their respective initial mid-plane 
densities, po, are given in Table I. The runs were started in hydrostatic 
equilibrium, using the density stratification profile of Ferriere (1998), 
but scaled by the density po. After of order 20-40 Myr of evolution, 
the system attains a convincing statistical equilibrium. The time se- 
quences shown in Figure 1 show the variations of the instantaneous 
SN II frequency after this initial period. 

Table I gives the average frequencies z^sn arising from our SN im- 
plementation. The value for the Average case, i^snii — 40kpc~ 2 Myr -1 
(which, naively extrapolated, gives 1/33 yr _1 for the whole Galaxy), 
was obtained after calibration of r c to ensure just such a suitable 
average rate. This value clearly cannot then be used in vindication 
of our model. The rates for the other two simulations subsequently 
arise with no further input, however, and the resultant variation is well 
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Table II. Two models with varying initial mid-plane magnetic 
field, B . Filling factors are defined as in Table I. The initial 
mid- plane density is po = 0.7 x 1CT 24 gcm~ 3 in both cases. 





Unit 


No field Strong field 


Initial mid-plane field, Bo 


fiG 


6 


Hot gas filling factor 




0.19 0.12 



fitted as u SN = 20kpc" 2 Myr -1 (p / 1(T 24 gem" 3 ) L7 (see Figure 1). 
We may consider our SN rate as a proxy for star formation rate (SFR), 
and the exponent above is consistent with observed scalings of the 
latter quantity with gas surface density S: SFR oc £ K , with k = 0.9— 
1.7 (Kennicutt, 1998). (For the density variations we impose, S scales 
directly with pq.) Whilst it is not clear that the two relations are directly 
comparable — the SFR should depend on column-integrated density, 
and our SN mechanism is more dependent on local volume density 
- this result remains of great interest, and potentially constitutes a 
meaningful test of the implementation of physics in our model. Further 
work on this correlation is planned. 

These simulations give volume filling factors of hot gas which are 
rather low (see Table I); plausible values for the Average case are of 
order 0.15-0.20. This is a facet of the model that is particularly sensitive 
to magnetic fields, however. The runs described above had an initial 
azimuthal field of the form B y = Bq cosh~ 2 (z/0.3kpc), with Bq = 6/iG. 
This Bq is a reasonable value for the total field in the solar vicinity, but 
that field should be split roughly equally between ordered and small- 
scale components. A uniform field is extremely effective in confining 
the expanding hot gas produced by SNe, as is evident from Figure 2, 
which shows horizontally averaged filling factors for two comparable 
runs with and without such a field. (We should comment that these 
filling factors may be rather sensitive to the precise definition of hot 
gas used, however. Also our use of closed boundaries in z may influ- 
ence the vertical structure of this quantity; future work should remove 
this latter limitation.) The difference in mean filling factor is almost a 
factor of 2, as can be seen in Table II. (A solution for an intermediate 
Bq, not continued for so long, gives consistent results.) Re-considering 
the values in Table I (obtained with the over-strong large-scale Bq) in 
this light, we believe that our model represents the evolution of such 
structures in the gas in a reasonable way; and we emphasise the clear 
importance of magnetic fields — neglected in most other simulations - 
in this respect. 
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Figure 2. Hot gas filling factor /, as a function of height and time, for runs with 
no magnetic field (left panel) and with Bq = 6 /iG (right panel) . Filling factors 
range from to 0.93; higher values are indicated by lighter shades. Both runs are 
for Interarm densities. 
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